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Abstract. The high order fluid model developed in the preceding paper is employed 
here to study the propagation of negative planar streamer fronts in pure nitrogen. The 
model consists of the balance equations for electron density, average electron velocity, 
average electron energy and average electron energy flux. These balance equations have 
been obtained as velocity moments of Boltzmann's equation and are here coupled to the 
Poisson equation for the space charge electric field. Here the results of simulations with 
the high order model, with a PIC/MC (Particle in cell/Monte Carlo) model and with 
the first order fluid model based on the hydrodynamic drift-diffusion approximation 
are presented and compared. The comparison with the MC model clearly validates 
our high order fluid model, thus supporting its correct theoretical derivation and 
numerical implementation. The results of the first order fluid model with local field 
approximation, as usually used for streamer discharges, show considerable deviations. 
Furthermore, we study the inaccuracies of simulation results caused by an inconsistent 
implementation of transport data into our high order fluid model. We also demonstrate 
the importance of the energy flux term in the high order model by comparing with 
results where this term is neglected. Finally, results with an approximation for the 
high order tensor in the energy flux equation is found to agree well with the PIC/MC 
results for reduced electric fields up to 1000 Townsend, as considered in this work. 
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1. Introduction 

Streamer discharges occur in nature and technology, predominantly in pulsed high 
voltage discharges; at their rapidly propagating tips electric fields well above the 
break-down value are maintained that create high electron energies and reaction rates. 
Streamers occur in lightning and sprites [H El [3l H] as well as in industrial applications 
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such as lighting El E], treatment of polluted gases and water [S], disinfection [H], 
plasma jets and bullets [101 IHl 1121 HSl E] and plasma assisted combustion [T6] . 
Further optimization and understanding of such applications depends on an accurate 
knowledge of the electron dynamics during streamer development. 

A Monte Carlo technique is able to track the electron distribution in phase space, 
yielding both electron density profiles and electron energy distribution. However, a full 
particle model long has been computationally too demanding or too inaccurate due to 
the used super-particle techniques, and even now the parameter range of simulations is 
limited. For recent progress in particle and hybrid models we refer to [HI [191 1201 [HI [22] . 
Due to the computational costs and limitations of particle models, up to today mainly 
fluid approximations have been used to model the structure and evolution of streamer 
discharges in two or three spatial dimensions. These fluid models were restricted 
almost exclusively to the hydrodynamic reaction drift diffusion approximation combined 
with the local fleld approximation, though the short-comings of this model have been 
acknowledged and documented [I8l[l9l[2a[23l[2a[25l[26l[271[28l[29l[30^ They arise 
from the fact that within the streamer front the electron density develops large spatial 
gradients and a complex interaction with the fleld. This in turn limits the use of the 
local fleld approximation as the electron energy does not immediately relax to the value 
determined by the local fleld, but depends upon the electric fleld in a wider spatial 
range. 

The aim of the preceding and the present paper is to develop and test a better 
fluid approximation for streamer discharges in gaseous media. In the preceding 
paper |3T] we have derived a high order fluid model, whose predictions for streamer 
dynamics will be tested in the present paper. The high order fluid model contains 
equations for the electron density, for the average electron velocity, for the average 
electron energy and for the average electron energy flux. It was obtained from velocity 
moments of the Boltzmann equation and closed in the local mean energy approximation. 
Momentum transfer theory [32l [331 [M] was used to evaluate the coUisional terms in the 
balance equations, and particular emphasis was placed upon the correct representation 
of momentum and energy transfer in energy-dependent non-conservative collisions. 
The system was truncated at the level of the energy flux balance and simplifying 
approximations were introduced in the momentum and energy flux balance equations 
in order to close the system. In particular, in the momentum balance equation it 
was assumed that the distribution of velocities is isotropic assuming isotropy of the 
temperature and pressure tensors. Thus, the pressure tensor was reduced to a scalar 
kinetic pressure. Furthermore, the same approximation was used in the energy flux 
balance equation. If one assumes an isotropic distribution of velocities then the 
contribution of higher terms including the high order heat flux tensor and the high 
order pressure tensor is relatively small and could be neglected. However, the energy 
flux transported by the convective particle motion should be treated and implemented 
carefully. The high order tensor in the energy flux balance equation is expressed in terms 
of lower moments using the simplifying assumption that the pressure tensor is isotropic. 
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Thus with such a treatment of the random motion and the high order terms associated 
with the energy flux of the drift motion, we have obtained a complete and closed system 
of fluid equations sufficient to determine all macroscopic streamer properties. 

In the present paper the system of fluid equations derived in the preceding paper 
is solved numerically for a negative planar streamer ionization front and compared with 
results of a PIC/MC particle model and of the first order fiuid model. The comparison 
with the particle model confirms that our high order fiuid model approximates the 
particle behavior very well. The comparison with the commonly used first order fiuid 
model of reaction drift diffusion type shows the short-comings of this classical model. 
This concerns the front velocities, the ionization levels behind the front as well as the 
electron energy distribution in the high field region ahead of the front and in the 
ionized low field region in the streamer interior. Actually the inaccurate ionization 
level in the streamer interior in the first order fiuid model is a direct consequence 
of the inaccuracies of the electron energy distribution [HI [19]. Inaccurate electron 
densities and energies also have direct consequences for the gas chemistry in the streamer 
interior |35l |36|, [371 [SHI ISH]- In addition to the accuracy and comprehensiveness of our 
high order fiuid model, its firm theoretical foundation is certainly an additional favoring 
factor. 

We begin by presenting in section [2] the differential equations governing the 
electron density, the average electron velocity, the average electron energy and the 
average electron energy flux. Then we discuss the numerical algorithm used for the 
solution of the differential equations. The speciflc elements of this discussion include a 
brief mathematical reminder on hyperbolic systems, the implementation of initial and 
boundary conditions and the discretization in space and time. Our numerical analysis 
is performed in a ID model, because in the present paper we are not interested in 
a streamer morphology and related phenomena where a full multidimensional model is 
required. We use a PIC/MC method as an alternative technique to verify our high order 
fluid model. Following the recent work of Li et al. [HI [191 120] , in section [3] we give a 
brief description of this method with particular emphasis upon the construction of planar 
fronts in the particle model. The high order fluid model is then applied in section[4|where 
examples of planar fronts in pure nitrogen are presented. The effects of photoionization 
are not considered and present theory can be applied to negative fronts in discharges 
where photoionization is very weak (high-purity nitrogen is a good example). Results 
demonstrating the validity of a high order fluid model compared with the PIC/MC are 
shown. Following the preceding paper, we demonstrate the importance of a consistent 
implementation of transport data in fluid models, but now instead of using the flrst 
order model as in the preceding paper, we employ our high order model. Along similar 
lines, the accuracy of the two term approximation for solving Boltzmanns equation in 
the context of high order fluid studies of the streamer discharges is examined. We pay 
particular attention to the role of the electron energy flux. For illustrative purposes, 
calculated results for the planar fronts based on a high order model with energy flux 
and those calculated without energy flux, are compared over a range of electric flelds in 
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order to verify the physical arguments associated with the solution regimes and closure 
assumptions outlined in the preceding paper. Next, we examine the closure assumption 
associated with the explicit influence of the high order tensor appearing in the energy 
flux equation on the streamer dynamics. We conclude that the explicit contribution of 
the high order tensor in the energy flux equation is negligible and it makes a notable 
difference only in the velocity of the planar fronts. In summary, in this paper we 
present the numerical solution of a high order fluid model, test and verify it on the more 
microscopic PIC/MC model, and discuss how inherent streamer properties deviate when 
simpler fluid models are used. 



2. Model description and numerical solution 

2.1. The high order fluid model 

In our preceding paper [31] we have developed a high order fluid model for streamer 
discharges. The balance equations were obtained as velocity moments of the Boltzmann 
equation while the coUisional terms were evaluated using momentum transfer theory. 
We have truncated the system of fluid equations at the level of energy flux balance by 
reducing the pressure tensor to a scalar kinetic pressure and by neglecting the high order 
terms associated with the flux of thermal motion. The energy flux of the drift motion, 
however, is included. For more details of the derivation we refer to [3.1j. The derived 
model consists of a set of differential equations for the electron density n, for the average 
electron velocity v, for the average electron energy e and for the average electron energy 
flux ^: 

— + V-nv = -n{uA - i^i) , (1) 
^{r^v) + A v(ne) - n^E = -nv (u^ + + ^e^) , (2) 

— (ne) + V ■ «) - eE ■ (nv) = -nu, (e - -kTo 

-nY^{V,- V:)e, - neV^ - nj^ ^ Ae« , (3) 

i i 

+ - ^neeE = -u^n^ , (4) 

where E is the electric field, m and e are electron mass and charge, Tq is gas temperature 
and k is the Boltzmann constant. The collision frequencies for momentum u^n and energy 
z/g transfer are given by equations (49) and (50) of the preceeding paper while and 
ua are the ionization and attachment rate coefficients. Ui and z7/ are inelastic and 
superelastic collision frequencies, respectively for inelastic channel i. /3 is a parameter 
introduced to approximate the high order tensors in the energy flux equation in terms 
of lower moments (see equation (55) of the preceeding paper). 
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Taking into account the rapid growth and propagation of the streamers with orders 
of the electron drift velocity and the relatively slow drift velocity of ions, the ions are 
approximated as immobile, as usually done when streamer ionization fronts are analyzed. 
Therefore, the charge densities due to the positive and negative ions change only due to 
ionization and attachment: 



-n[UA - yj) 



(5) 



In order to account for the space charge effects, the system (P)-© is coupled with the 
Poisson equation for the potential which then enables the calculation of the electric 
field E: 



— {ni. 



n) 



E 



(6) 



where eo is the dielectric constant. 

2.2. ID hyperbolic system of balance laws 

In the present paper we simulate the propagation of negative streamer fronts in pure 
N2 at atmospheric pressure and at an ambient temperature of 298 K in ID. As N2 is 
a non-attaching gas, electron attachment does not appear in the source terms of our 
system ([I])-©. The effects of superelastic collisions are not taken into account and 
we consider only single ionization with ionization energy ej. The electron dynamics of 
equations in ID is given by the following nonlinear system of balance laws: 



du ^ , .du 



F(u) 



where the primitive variables are 

u = (n, nv, ne, )'^, 
the matrix Aiu) is 
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and the source term is 
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The equation for the ions (I5j) is without electron attachment 



dt 
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while the Poisson equation (jS]) in ID has the following simple form: 

^ = —{riion - n) . (12) 
ox 60 

Hence, the streamer dynamics in ID for a non-attaching gas is described by the system 
of equations fl71)- flT^ . 

Now we can apply the theory for hyperbolic systems of balance laws [HI HJl US] to 
our system (I71)-(l9]). The matrix A{u) has four eigenvalues: 

Ai,2 = ±7^ /3 - ^m^) ^^3^ 
A3,4 = ±7^/3 + ' 
where 7 = -y/^- AH eigenvalues are real and distinct when 

/3 = or 13 > I. (14) 

This means that the system ([7]) is hyperbolic if the condition f|T^ holds. Although 
the eigenvalues ( |T3l) are simple, the corresponding right and left eigenvectors have a 
rather complicated form, which makes it impossible to work with the vector of Riemann 
invariants. 



2.3. Initial and boundary conditions 

The electric field E = Ee^ (where is the unit vector in the x direction) drives 
the dynamics. We take as a positive value; therefore electrons drift to the left, and 
negative streamer ionization fronts move to the left as well. To create steady propagation 
conditions for the negative front, the electric field on the left boundary x = is fixed to 
the time independent value Eq 

E{0,t) = Eo>0. (15) 

The electric field for x > is calculated by integrating equation (fT2il numerically over 
X, with (fT5|) as a boundary condition. The right boundary is located at x = L; in all 
calculations we set the system length L to 1.2 mm. 

All simulations are started with the same initial Gaussian type distribution for 
electrons and ions 



n{x)\t=o = riiexp 



{x - Xpf 
a2 



(16) 



where we have chosen rij = 2 x 10^^ m^^, xq = 0.8 mm and a = 0.029 mm. The initial 
conditions for the average electron velocity, average electron energy and average electron 
energy flux are taken to be spatially homogeneous. The actual values of these quantities 
are calculated using a multi term solution of Boltzmann's equation, as already discussed 
in the preceding paper. 

As discussed in the appendix, we use homogeneous Neumann boundary conditions 
for all components of u at both ends of the system, so that all electron quantities may 
flow out of the system. Only for the electron density at x = L we employ a homogeneous 
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Dirichlet boundary condition to prevent electrons from diffusing out. However, all 
calculations end before the ionization reaches the boundary, that means, before the 
boundary conditions start to become relevant. 

Further details on the numerical discretization of the system can be found in the 
appendix. 



3. MC particle model and first order fluid model 

In the next section, we will compare the simulation results of our high order model with 
those of the PIC/MC particle model and with the first order fluid model. Here these 
models are briefly described. 



3.1. The MC particle model 

A Monte Carlo model following the motion of individual electrons contains the full 
physics that is to be approximated by a fluid model. The successful comparison of a 
fluid model with such a particle model validates the fluid model; of course, the same 
cross-sections have to be used in both models. 

The construction of a planar front is straight forward in fluid models, as the 
spatial derivatives are simply evaluated in one direction. However, in the particle model 
electrons move in all three spatial dimensions and hence, the three-dimensional setting 
has to be restricted as in previous work [18]. An essentially one- dimensional setting is 
achieved by considering only a small transversal area T of the front and by imposing 
periodic boundary conditions at the lateral boundaries. The electric field is calculated 
only in the forward direction x through 

E,{x,t) = E,{xo,t) + dx' ^- "--^ )_^^y^^^ (^7^ 

Jxo JT T eo 

where is the electric field in the x-direction, and xq is an arbitrary position. Therefore, 

fluctuations of the transversal field due to density fluctuations in the transversal 

direction are not included. The density fluctuations projected onto the forward direction 

depend on the transversal area T over which the averages are taken. 



3.2. The first order fluid model 

The first order fluid model is the multiply used reaction drift diffusion model that also 
was called the "classical fluid model" in our previous work [HI [22]; in these previous 
papers, it was already shown that the model can only serve as a first approximation, 
but cannot reproduce the results of a particle simulation quantitatively. 

The steps of derivation of this model were discussed and criticized in our 
previous paper [31]. The model is based on the hydrodynamic reaction drift diffusion 
approximation 

dn 9 / _ _(9n\ , 
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where /z and D are electron mobility and diffusion coefficient, respectively. Transport 
coefficients used as an input in this model are functions of the local electric field as 
discussed in the preceding paper. In order to account for the space charge effects, 
equation (fT8|) is coupled to the equations for ion density (fTT]) and for the electric field 
( IT2I) . For the purpose of comparison, we use the same numerical schemes to discretize 
equation f|T8|) and the same initial and boundary conditions as introduced above for the 
high order model. 

4. Results and discussion 

In this section we present and compare the simulations results obtained with our high 
order fluid model, with the PIC/MC method and with the first order fluid model. We 
consider the configuration described in Section 12.31 The transport coefficients required 
as an input in both the first and the high order fluid model are given in our preceding 
paper [31] . They were calculated from cross sections for electron scattering in N2 through 
a multi term solution of the Boltzmann equation. Details of the calculation together 
with the prescription how to use the data in modeling are given in the same article. It 
is important to emphasize that the same coUisional cross sections for electrons in N2 are 
used in all models presented in this paper. 

4.1. Simulation of planar fronts: overall comparison between the models 

Figured] shows the temporal evolution of the electron and ion densities and of the electric 
field in N2, when the electric field ahead of the front is 590 Td (or equivalently 145 kV/ cm 
at atmospheric pressure and at a temperature of 298 K). The initially Gaussian electron 
density shows the behavior observed many times in past: first, it grows due to the 
ionization processes; then charge separation occurs in the electric field due to the drift 
of oppositely charged particles in opposite directions, and the initially homogeneous 
electric field is distorted; and finally, when the field in the ionized region becomes more 
and more screened, the ionization stops in the ionized region and the typical ionization 
front profiles of electron and ion densities and of the electric field are established. The 
initial ionization avalanche is then said to have developed into a streamer. 

We note that the results of the high order fluid model and of the PIC/MC simulation 
agree very well. The streamer velocity is almost the same, as well as the electric field 
and the electron/ion density. Conversely, when the first order fluid model is applied, the 
streamer velocity is lower, with a lower electron density in the streamer head and in the 
streamer channel. In figure |2] we show how the streamer velocity depends on the applied 
reduced electric field E/uq. In order to calculate the streamer velocity we have followed 
the evolution of a certain level (2nj) of the electron density at the streamer front. We see 
that for increasing E/uq the differences of the streamer velocities between the first order 
and the high order fluid model increase slightly. On the other hand, the velocities of the 
high order fluid model and the PIC/MC agree very well. Similar trends are observed for 
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Figure 1. Temporal evolution of the electron density (a), ion density (b) and electric 
field (c) in a planar front in N2. Shown are the spatial profiles obtained with three 
different models as indicated on the graph. Flux transport data are used as an input 
in fluid models. The externally applied reduced electric field is 590 Td. 



the electron density behind the ionization front. While the first order model generally 
underestimates this electron density, the high order model approximates the PIC/MC 
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results much better. This agreement holds in the full range of fields E/uq explored, 
while the differences between first order model and MC model increase with increasing 
E/no. 

In some aspects the results presented in figures [1] and [2] are consistent with 
observations by Li et al. [iHl [19] and by Kanzari et al. [28]. In particular, Li et 
al. [T9| [22] already have found considerable differences in velocities and ionization 
densities between MC particle model and classical fluid model; they have proposed an 
extension of the first order fluid model in [19], based on a phenomenological gradient 
expansion suggested in [531 El] to get a better agreement of the fluid model with the 
results of the particle model. The extended fluid model of Li and co-workers uses a 
density gradient expansion of the source term to approximate the spatial non-locality 
of the ionization processes at the streamer front. However, this model overestimates 
the particle results for the electron density in the streamer channel when the field is 
below 125 kV/cm (for standard temperature and pressure of the background gas) by 
up to 4 %, while for an electric field of up to 200 kV/cm, the extended fluid model 
underestimates the particle results by up to 9 % [19] . On the other hand, the differences 
in the electron density behind the ionization front between our systematically derived 
high order fluid model and the PIC/MC model do not exceed 5% for electric fields of up 
to 245 kV/cm as considered in this work, thus demonstrating that we correctly derived 
and implemented our high order fluid model. 



— o— high order bull< 
— *— first order bulk 
-□-PIC/MC 




0.2-^ — , — ■ — , — ■ — I — ■ — , — ■ — I — ■ — I — ■ — , — ■ — I — ■ — I 
300 400 500 600 700 800 900 1000 1100 

E/n, [Id] 

Figure 2. Velocities of planar fronts as a function of the electric field obtained with 
the first and the high order fiuid model and with the PIC/MC method; different sets 
of input data are used as indicated on the graph. The electron flux drift velocity as a 
function of E/tlq is also included. We remark that the negative sign of all velocities is 
removed here. 
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4.2. Electron energy profiles in different front regions within the different models 

We now focus on the characteristic differences in the profiles of the electron density. 
Figure [3] shows the profiles of the average electron energy for all three models, as well as 
the profiles of the electric field and of the electron density at t = 0.7 ns for the first order 
model only. In the first order model, the average electron energy is assumed to adapt 
to the electric field instantaneously, therefore it was derived from the local electric field 
for the plot. However, the profile of the average energy obtained by the high order fiuid 
and PIC/MC models is more complex. 




Figure 3. The average electron energy in three different models at time 0.7 ns. Profiles 
of the electric field and the electron number density from the first order fluid model 
are also included. Calculations are performed for E/no of 590 Td. 

We distinguish three different spatial regimes, the streamer head where most 
ionization occurs and where electron density and electric field have large gradients, 
the region ahead of the front where the electric field is high and the electron density 
vanishes, and the streamer interior where the electron density is finite and the electric 
field is low or vanishing. 

First, figure [3] shows that in the streamer head the average electron energy given 
by the high order fluid and PIC/MC methods are higher than in the first order fiuid 
model. Since the ionization rate in this energy range is an increasing function of the 
electron energy, it is clear that the higher the average electron energy, the higher is the 
ionization rate. These observations are consistent with those made by Li et al. p!8| 119]. 

Second, in the region ahead of the front the average energy in the high order model 
has a slope which reflects the so-called non-local effects. In this region, the PIC/MC 
method is not efficient as there are not enough particles to sample the spatially resolved 
average electron energy correctly. It must be emphasized that the spatial variation of the 
average electron energy is present even in the avalanche phase where an initial spatially 
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homogeneous electric field is not distorted due to space charge effects. According to 
arguments given in fL8\, the leading edge of an ionization front has essentially the 
same dynamics and therefore also the same energy slope as an electron swarm. In 
the past, many swarm oriented studies have been performed to explore the effects 
of spatial variation of the average energy through the swarm and many phenomena 
have been discovered [SU HU HHl |16]. Almost certainly, one of the most striking 
phenomena is the difference between the flux and the bulk transport data which follows 
directly from spatially dependent nonconservative collisions (ionization and/or electron 
attachment) resulting from a spatial variation of average electron energies within the 
swarm [Ml El HU HE] • 

Third, figure [3] shows strong disagreement between the average energy in the first 
order model and the average energies predicted by high order fluid and PIC/MC methods 
in the streamer interior. Within the first order fluid model the average electron energy 
attains approximately the thermal value. This follows directly from the local field 
approximation and from our Boltzmann equation calculation of the mean energy in 
the limit of vanishing electric field. It should be emphasized that thermal effects of 
background molecules are included in our Boltzmann equation analysis. In contrast to 
the first order fluid model, the average electron energy given by high order fluid and 
PIC/MC methods are significantly higher. This is another manifestation of the non- 
local effects, but this time, these phenomena take place in the streamer channel. As 
the velocity of the front is higher than the electron drift velocity even in the streamer 
head region, electrons slowly relax to lower energies in the lower or vanishing field in the 
streamer interior. This causes the spatial decay of average electron energy in high order 
fluid and PIC/MC model, that can be observed in figure [S This decrease of the average 
electron energy follows directly from the energy dependence of the collision frequency 
for energy relaxation shown in figure 2 of the preceding paper [3T] . 

Another interesting feature is the disagreement between the high order and the 
PIC/MC results for the average electron energy in the streamer channel. We believe that 
this disagreement arises from the limitations associated mainly with the high order fluid 
model. Although our PIC/MC method does not treat thermal effects of the background 
molecules, electron energies are too high for this to be significant. Most probably this 
disagreement arises from an increased anisotropy of the distribution function which in 
turn limits the adequacy of the closure assumptions associated with the pressure tensor. 
Indeed, the average energy in the streamer channel varies between 0.5 and 1 eV and 
exactly in the same energy region the cross sections for vibrational excitation grow 
rapidly in magnitude (a few orders of magnitude in a very narrow energy range) while 
the cross section for momentum transfer in elastic collisions varies relatively slowly with 
the electron energy. This will induce strong anisotropy of the distribution function, and 
it is clear that the assumption of the isotropy of the temperature tensor is certainly 
problematic under these conditions. These physical arguments are verified on figure 
m where the ratio between transverse and longitudinal components of the temperature 
tensor for electrons in N2 is shown. It shows a deep minimum in the profile between 0.5 
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Figure 4. The ratio between transverse and longitudinal components of the 
temperature tensor as a function of electron energy in N2 based on the multi term 
approach for the Boltzmann equation. 



and 1 eV which is a clear sign of an increasing anisotropy of the distribution function in 
velocity space. It should be emphasized that for the range of electric fields considered 
in this work, the different average energies in the streamer channel observed in high 
order fluid or PIC/MC model cannot induce significant changes in the electron density. 
Another physical argument which can be used to address the differences between the 
average electron energies in the streamer channel is associated with the accuracy of 
the PIC/MC method. In the energy range around 1 eV one must carefully follow the 
electrons in a PIC/MC simulation. If the time step between two successive collisions 
is too big then the energy losses due to rapidly increasing cross sections for vibrational 
excitation are not going to be included accurately. This is of less importance for the 
coUisional processes whose cross sections do not vary so rapidly with the electron energy. 
One way to overcome this problem is to reduce the time step in a PIC/MC method but 
the penalty is a significant increase in the computation time. It is clear that additional 
testing and calculations are required to resolve this issue. 

Figure [5] shows the average electron velocity through the streamer front. As for 
the average electron energy, in the first order model the average velocity of electrons is 
obtained from the profile of the electric field assuming the local field approximation. 
We see that the average velocity of electrons follows changes in the electric field 
instantaneously: in the outer streamer region where the electric field is constant the 
average velocity has also a constant value while in the streamer channel the average 
velocity is essentially zero. When considering the profiles of the average velocity derived 
by the high order model, one may observe the characteristic slope of this quantity in 
the outer region of the streamer while in the streamer channel it vanishes. This stands 
in contrast to the profile of the average energy. The relaxation of the average velocity is 
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Figure 5. Profiles of tlie average electron velocity in the first and the high order fluid 
models and of the PIC/MC model at time 0.7 ns. Profiles of the electric field and of 
the electron number density from the first order model are also included. Calculations 
are performed for E/uq of 590 Td. We remark that all velocities are negative, and 
their absolute value has been plotted. 

determined by the momentum transfer collision frequency which is a factor of 2.5 x 10^ 
larger than the collision frequency for energy transfer. This means that the relaxation 
of the average velocity is much faster and hence in contrast to the average electron 
energy this quantity relaxes very rapidly. Figure |5] shows a good agreement between 
results derived by fluid models and by the PIC/MC method. As in the case of the 
average electron energy, the PIC/MC method is not suitable for the determination of 
the average velocity of electrons in the region ahead of the front because in this region 
only a few electrons exist. 

4.3. On the use of transport data in the high order fluid model 

Another important issue in streamer modeling concerns the choice of transport data 
in fluid models. The origin of, and the difference between bulk and flux transport 
coefficients in the context of streamer studies have been recently discussed in p^[20llTf] . 
In this section we consider the implications of the transport coefficient duality in the 
context of our high order fluid model for streamers. We also discuss to what extent 
the differences between transport coefficients obtained by the two term approximation 
for solving Boltzmann's equation and those obtained by a multi term theory affect the 
accuracy of streamer models. 

Figure [6] shows the temporal evolution of the electron and ion densities and of the 
electric field in a streamer front when the reduced electric field ahead of it is 590 Td. 
The calculations are performed for three different sets of input data as indicated in the 
figures. The inadequacy of using bulk data in the high order model is clearly evident 
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Figure 6. Temporal evolution of (a) the electron density, (b) the ion density and (c) 
the electric field in a planar negative ionization front in N2 with a reduced electric field 
E/no of 590 Td ahead of the front. Displayed are the spatial profiles obtained with 
three different sets of input data as indicated in the graph. 



High order fluid model for streamer discharges. II. Investigation of planar fronts. 



16 



flux bulk bolsig+ 




0.0000 0.0002 0.0004 0.0006 0.0008 




c) 




0.0008 



Figure 7. The same simulations and plots as in the previous figure, but now for (a) 
the average electron velocity v, (b) the average electron energy e, and (c) the average 
electron energy flux ^. Again, velocity v and energy flux ^ here are negative quantities, 
and their absolute value is plotted. 
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(having in mind the very good agreement between the profiles obtained with the flux 
data and with the PIC/MC method, as shown in the previous section): the streamer 
velocity is higher and the electron/ion density in the streamer head and in the streamer 
channel is higher than with the flux transport data. For increasing electric fleld the 
discrepancies in the streamer velocity increase further, as illustrated in flgure O Similar 
trends have been observed in the preceding paper |31] where planar fronts were modeled 
with the first order fluid model. In both models, streamers with bulk data are faster and 
the electron/ion density is higher in both the streamer head and the streamer channel. 
This follows directly from the fact that in N2 the bulk mobility of the electrons is larger 
than the flux mobility. On the other hand, as already shown in section 14.11 we observe 
an excellent agreement between our high order proflles obtained with the flux data and 
those predicted by a PIC/MC method. 

Figure [7] displays the temporal evolution of the average electron velocity, of the 
average electron energy and of the averaged electron energy flux under the influence of 
electric fleld of 590 Td for different times as indicated on the graph. The inadequacy of 
using bulk data is again evident. In the early stage of the streamer development (when 
the electric fleld is not entirely screened), the average electron energy calculated using 
the bulk data is higher in all relevant streamer regions. When the transition process 
from an avalanche to a streamer is flnished, we see a clear discrepancy between the 
average energies in the region ahead of the streamer front. As simulated planar fronts 
with bulk data are faster, the corresponding proflles of the average energy are shifted 
forward. However, for a fully developed streamer, there are no signiflcant deviations 
between the average energies in the streamer channel. In this region the electric fleld is 
screened, and the average energies are slowly thermalizing although not fully relaxed as 
shown in previous section. The average electron velocity behaves exactly in the same 
manner as the average electron energy when comparing results with flux or bulk data. 
The only difference is associated with the fact that the average velocity of electrons 
is almost fully relaxed in the streamer channel regardless of the type of data used in 
calculations. The average electron energy flux shows similar behavior. 

In earlier streamer investigations it has not been generally investigated to what 
extent the two term Boltzmann equation results for various transport coefficients affect 
the accuracy of the model predictions. The limitations of the two term approximation 
for solving Boltzmann's equation have been illustrated many times in past in the context 
of swarm studies [M l H91I501I5T] . Here the question arises of whether a similar conclusion 
can be drawn for streamers taking into account that the streamer development is a non- 
linear, non-stationary and non-hydrodynamic problem where space charge effects are 
important and where the electron energy varies from thermal values in the streamer 
channel up to a few tens of eV in the outer region of the streamer for the fields 
considered here. It is clear that the anisotropy of the velocity distribution function 
may be considerable in various streamer regions and hence we will compare profiles for 
different streamer properties using the two term and multi term results for transport 
data as an input in modeling. 
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Figure E] (a) shows that the ionization level behind the streamer front calculated 
using the BOLSIG+ transport data is lower while the streamer velocity is the same. 
Surprisingly, the profiles of the electric field, the average electron energy and the 
average electron velocity obtained with the multi term flux data or with the BOLSIG+ 
data agree very well. This is a clear indication that internal errors of the two term 
approximation associated with the ionization rate can be directly used to understand 
the differences in the profiles of the electron and ion densities in the streamer channel 
(see figure 1 (c) in the preceding paper). However, this conclusion cannot be generalized 
to other gases or other conditions than those studied here. As pointed out recently by 
White and co-workers [311 EO] one may never be sure about the accuracy of the two term 
approximation as various transport properties show different sensitivity with respect to 
this approximation in different energy ranges. For example, the inadequacy of the two 
term approximation has been recently demonstrated for CF4 [19] and CH4 [50]. In 
the preceding paper, it was shown that the velocity of a streamer in the first order 
fluid model is significantly affected by choosing either the two term or the multi term 
solution. Therefore, when high precision is required the best option is to use a multi 
term approach and/or a Monte Carlo simulation technique to calculate the input data 
for fluid models regardless of their order. 

4.4- Simulation of planar fronts with and without the energy flux term 




0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 



X r m 1 

Figure 8. Electron density, ion density and electric field in N2 in an external electric 
field of 350 Td. The two lines indicate solutions of the high order fluid model with 
and without the energy flux term at the same instance of time. The initial condition 
is the same, and the flux transport data are used. 

In this section, we compare the results of the high order fluid model with and 
without the energy flux term Q. In the model without the energy flux term the energy 
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Figure 9. The same instant and plot as in Fig. [8l but now average electron energy, 
negative average electron velocity and negative average electron energy flux are plotted. 

flux ^ in the energy balance equation is set to zero, and therefore the energy flux 
balance equation (jl]) does not need to be considered. This approach was aheady used 
in the case of streamer corona discharge dynamics [281 129] and for analysis of ionization 
wave dynamics in low-temperature plasma jets [13]. Figures El and [9] show the streamer 
properties for a reduced electric field E/uq of 350 Td at time t = 1.5 ns. First, we see 
that the high order fluid model with the energy flux term gives higher electron and ion 
densities within the streamer channel than without the energy flux. For increasing E/rio 
these differences increase further, as shown in figure [TOlfa). The electric field exhibits the 
typical streamer behavior; the observed differences between the high order fluid models 
with and without the energy flux term follow from the time delay needed for the space 
charge to become high enough and distort the externally applied spatially homogeneous 
electric field in these two models. This property is even more evident for the streamer 
velocity shown in figure [TOT b). Figures El and [9] show clearly that in our simulation 
conditions the streamer front propagates faster in the high order fluid model with the 
energy flux term. 

Furthermore, the average electron energy and the average electron velocity ahead of 
the streamer depends clearly on the fact whether the electron energy flux is included or 
not. The average electron energy has a steep gradient between the outer streamer region 
and the streamer channel. In this spatially narrow region the energy decreases from the 
range between 10 and 20 eV to essentially thermal values in the streamer channel. It 
is clear that the correct treatment of the energy flux in the high order fluid model is 
critical for an accurate energy transport between the different streamer regions. In view 
of the close agreement between the present results with the energy flux term and those 
predicted by a PIC/MC method, established in the previous sections, we must conclude 
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Figure 10. (a) The ionization level behind the planar front, and (b) the absolute 
value of the front velocity as a function of the reduced electric field EjriQ ahead of the 
front. Results are shown for the high order fluid model with or without the energy flux 
term, as well as for the PIC/MC model. Flux transport data are used as an input in 
the fluid models. 



that the high order fluid model without the energy flux term is qualitatively right but 
quantitatively wrong under the conditions simulated. 

Simulation of planar fronts with and without the high order tensor 

In this section we investigate the influence of the high order tensor in the energy 
flux equation. In order to truncate the system after the energy flux equation, in the 
preceding paper this tensor was expressed by a product of two lower order tensors times 
a parameter (3. In all previous sections, we used /? = 1 which is a natural guess. Here 
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we investigate how the simulation results depend on /3. 
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Figure 11. Profiles of (a) electron density, (b) electric field, (c) average electron 
energy, and (d) average electron energy flux tor (3 — I, (3 — 2 and /? = 3. The results 
are shown for E/uq of 350 Td at time 1.5 ns. 

Figure [11] shows the profiles of the electron density, the electric field, the average 
electron energy and the average electron flux for a reduced electric field E/hq of 350 Td 
at time 1.5 ns. One can immediately see that increasing the parameter P increases the 
velocity of the streamer very slightly. All other quantities remains essentially unaltered 
in the streamer channel indicating their weak sensitivity with respect to /3. Only in the 
limit of the highest electric field considered here, the streamer properties in the channel 
respond slightly to the variation of the parameter /3. Ahead of the streamer front, 
however, we observe variations in the streamer properties even for the reduced electric 
field of 350 Td. As expected, the most sensitive quantity with respect to variations in 
P is the average electron energy flux. 

In section 12.21 we have derived that our high order model is well posed if /3 = or 
/3 > 1. We have solved our system for a range of values for /3 and reduced electric fields 
E/uq and it was found that (3 = 1 ensures the best agreement between profiles obtained 
by the high order fluid model and by the PIC/MC method. This validates the closure 
assumption associated with the contribution of high order tensors in the energy flux 
equation, as discussed in the preceding paper. Furthermore, if we analyze the profiles of 
various streamer properties for /3 = 0, the simulation results agree very well with those 
derived by the PIC/MC method. The importance of this is twofold. First, the influence 
of the high order tensors in the energy flux equation is less than initially expected. 
Second, taking /3 = is a very good approximation, and it will significantly reduce the 
computation costs for the differential equations in 3D. 
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5. Conclusion 

In the previous paper we have derived the equations and the transport and reaction 
parameters of a high order fluid model for streamer discharges. In the present paper 
we have briefly presented an accurate and efficient way to numerically solve the high 
order fluid model. These numerical schemes are then used to study the propagation 
of negative streamer fronts in N2. Then we compared the high order fluid model and 
the flrst order fluid model with a PIC/MC model. While the first order model is very 
commonly used for streamer simulations, it can be clearly seen that the high order model 
is a much better approximation of the full particle dynamics. This concerns not only 
quantitative, but also qualitative differences. In particular we observe: 

(a) The non-local effects in the profile of the average electron energy are present both 
in the PIC/MC-model and in the high order model, but are missing in the first order 
model due to the local field approximation. The slope of the average electron energy 
in the region ahead of the streamer is related to the spatial variation of the average 
electron energy in the avalanche phase of the streamer development. And the slope of 
the average electron energy in the streamer channel is striking and an inherent property 
of a streamer as well, as the electrons can not immediately relax to the low or vanishing 
field in the streamer interior. 

(b) Bulk and fiux transport coefficients for electrons can differ substantially in so-called 
non-conservative regions where ionization, attachment, detachment and recombination 
takes place. Here one needs to take care which coefficients to implement. Streamers 
in nitrogen with bulk data are faster and create a much higher electron density in the 
streamer channel. When the streamers in nitrogen are calculated with BOLSIG-I- data, 
their velocity is not affected. However, in the streamer channel and in the streamer head 
we have found differences between the electron densities in our high order fluid model 
when we inserted either BOLSIG-I- data or our flux data for the electron transport 
coefficients. 

(c) When we conducted streamer simulations with our high order model, either including 
or excluding the energy fiux, various streamer properties showed that the energy fiux 
term cannot be neglected. Streamers without energy flux term are slower and create a 
lower electron density in the streamer channel. In the region ahead of the streamer, the 
average electron energy and the average electron velocity also change when the energy 
flux is neglected. 

(d) The validity of the closure assumption associated with the treatment of high order 
terms in the energy flux equation is studied through the variation of the parameters 
used to parameterize the expression in which high order tensors are expressed in terms 
of lower order moments. We have pointed out that although the average electron energy 
and average electron energy flux in the outer region of the streamer are influenced by this 
parameter, the demand for a strict consideration of the high order tensor in the energy 
flux equation may be relaxed. This follows from the fact that it is almost impossible 
to notice the differences in the streamer velocity and electron density in the streamer 
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channel for different parameters used to parameterize the high order terms. 

Though in the present paper, we have tested our model on planar fronts, the theory 
of the high order model has been developed in the previous paper without restrictions 
for the full three-dimensional setting. An important aspect to be addressed in the future 
is our closure assumption for the pressure tensor that we assumed to be isotropic. This 
restriction will be removed in future work. 

Finally we would like to emphasize that though this study has concentrated on 
electron dynamics mainly in planar streamer fronts, the theory and the associated 
numerical solution of the system of differential equations are equally valid for ions up 
to the energy balance equation. Therefore we have derived a new modeling frame work 
that is widely applicable in reactive plasmas in full 3D. 



Acknowledgments 

It is a pleasure to acknowledge the helpful discussions with R.D White, W.H. 
Hundsdorfer, J. Teunissen and J. Bajars. AM and SD acknowledge support from 
STW-project 10751, SD acknowledges support from STW-project 10118, part of the 
Netherlands's Organization for Scientific Research (NWO). SD has also been supported 
by the MNTR, Serbia, under the contract number ON171037. 



Appendix 

In this appendix we briefly describe the discretization of our system of differential 
equations. The finite volume method is used to discretize the system fl71)- f|T2l) in 
space. The basic principle of the control volume method is to maintain the conservative 
properties of the system (JT]) over every volume element. We introduce control volumes 
or cells Vj as follows: 

Vj := [jAx, (j + l)Ax] , X, := (^j + Ax, j = 0, 1, M - 1, (A.l) 

where Ax = L/M is the spatial grid size while L is the length of the simulation domain. 

The numerical solution Uj{t) has to be interpreted as an approximation of the 
average value of n(x, t) over the control volume Vj at time t, i.e. 

Uj{t) = ^ j u{x,t)dx. (A.2) 

To approximate the spatial derivative in ((Tj) we use the second-order central 
difference discretization [52] . 

"^^l^' = ^— (u(x + Ax, t) - u(x - Ax, t)) + O(Ax^), (A.3) 
ox 2Ax 

which has already been used in previous first order fluid models and associated codes 
(see for example [IH])- Hence, the equation can be re- written in difference form: 



duj{t) 1 



dt 2Ax 



A{u,{t)){u,+r{t) - u,-i{t)) = F(u,{t)), (A.4) 
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where i = l,2, ....,M - 2. 

As discussed in the main text, homogeneous Neumann or Dirichlet boundary 
conditions are imposed: At the left boundary, i = 0, we impose a homogeneous Neumann 
boundary condition 

^^^ = or u^,it) = uo{t). (A.5) 
ox 

In this manner, we provide the value U-i(t) = |Ax, t) at i = — 1 which exceeds the 
computational domain. At the right boundary, i = M — 1, we impose a homogeneous 
Dirichlet boundary condition for the electron number density, while for the other 
components of u we impose homogeneous Neumann boundary conditions to obtain 
a value for UMif) = ""-(-^ + |Ax, t). 

To approximate the electric potential (|6]) we use the same strategy as for the 
primitive variables u, i.e. is averaged over control volumes. Hence, the electric field 
E{x,t) is discretized on the edges of the control volume, denoted by E{x,_i,t). But 
since we consider a ID case, the electric field follows directly from the charge densities 
by integrating equation (fT2|) along the x-direction: 

E{x, t) = E{0, t) + — ( [ (riionix', t) - n{x', iU Ax'. (A.6) 



^0 \J0 

We then can calculate E{xj,t) by: 

Eix,,t) = ^iEix^_i,t) + Eix^^,t)). (A.7) 

Now we are going to consider the ODE system flA.4l) . To approximate the time 
derivatives we use classical RK4 (Rounge-Kutta 4) time-integration scheme [S2], which 
is a fourth order method. This is an explicit method which always has a bounded 
stability domain. In our case the stability condition has the following form: 

max \XiAt/2Ax\ < C, (A.8) 

l<j<4 



or by taking into account /3 > y /3 + a/ /3{/3 — 1), when /3 > 1, we have: 



A^/2max£^ 

^2AxV 3m - ^ ^ 

This condition is called CFL stability criterion, where C depends on the particular 
time- integration method and space discretization. In our simulations, C is set to 0.1. 
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